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Abstract When highly viscous oil is produced at low temperatures, large pressure 
drops will significantly decrease production rate. One of possible solu¬ 
tions to this problem is heating of oil well by hot water recycling. We 
construct and analyze a mathematical model of oil-well heating com¬ 
posed of three linear parabolic PDE coupled with one Volterra integral 
equation. Further on we construct numerical method for the model and 
present some simulation results. 

Keywords: Oil well, integro-differential equation, Volterra integral equation 


Introduction 

An oil well producing at low temperatures may experience large pres¬ 
sure drops due to high viscosity of oil and wax forming. One way to 
avoid these pressure drops is heating of oil by hot water recycling. 

The tubing is surrounded by two annulus for water circulation. Hot 
water is injected into inner annulus and it flows out of the system through 
the outer annulus. The main technical concern is minimization of energy 
lost in the system while keeping oil temperature sufficiently high. 

Configuration just described will be called counter flow exchange. 
If the hot water is injected into outer annulus and leaves the system 
through inner annulus, then we talk about parallel heat flow exchange. 
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Figure 1. Counter flow heat exchange 


The outline of the paper is as follows. In the first section we present 
a simple one-dimensional mathematical model describing the heat ex¬ 
change in the system. We present only counter flow configuration since 
parallel flow configuration differs only in signs of water velocities. Solv¬ 
ability of a system of integro-differential equations describing the heat 
exchange is discussed in the second section. It is shown that the re¬ 
sult of Artola pQ can be applied. In final section we discuss numerical 
method for solution approximation and present some numerical results 
for counter flow and parallel flow configurations. 

A problem similar to this one was considered in engineering literature 
in [nj. 

1. Mathematical model 

Cross-sectional mean velocities of oil and water in inner and outer an¬ 
nulus will be denoted by Vo, Vi and Ve- They are assumed to be constant, 
and therefore the fluids have constant pressure drops. Furthermore, to 
simplify the model, we neglect friction and we take mass densities po 
(oil), pu) (water) to be constant. The heat is transferred between the 
tubing, inner and outer annulus and the formation according to New¬ 
ton’s law. 

With these simplifying assumptions and taking direction of z axis 
vertically downwards, we obtain the following three parabolic equations 
(see [2] for example): 
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for z G (0,-L) and t G (0,fmax)- The main variables are the temper¬ 
atures of oil, water in inner annulus, water in outer annulus and the 
temperature of the formation, denoted respectively by To, Tj, Tg and 
Tf. All coefficients are constant and they have the following meaning: 
ao = AoPoCo, ai = AiPujCw, ae = A^PwCw where Ao, Ai and A^ are 
cross-sectional areas and Cq, Cw are heat capacities. By bo, bg and bf 
are denoted heat transfer coefficients from Newton’s law, and by Do, 
Di and thermal conductivities, multiplied by cross-section areas. In 
counter flow exchange all three fluid velocities are positive. From mass 
conservation it follows aiVi = a^Ve- 

Heat flow in surrounding formation is assumed to be radial with re¬ 
spect to the tubing and to have constant (geothermal) gradient in ver¬ 
tical direction. We denote by Tz{z) geothermal temperature and by 
Ts{r, z, t) the temperature in the soil. Formation temperature 7/ is then 
given by Tf{z,t) = Ts{rf,z,t), where rf is formation radius. The tem¬ 
perature Tg is the solution of the heat equation with initial temperature 
Tz, temperature at infinity equal to T^, and prescribed heat flux qj at 
r = rf. In the other hand, qp is given by Newton’s law 


Qf = bf{Te-Tf). 


(4) 


Then, by applying Duhamel’s principle we can represent formation tem¬ 
perature by the formula 

/■* d 

Tf{z,t) =Tz{z) + J p{t-T) — qf{z,T)dT, (5) 

where p{t) = P{rf,t — T^j^i^kf {kp is thermal conductivity of the soil) 
and P{r, t) is the solution of the problem 


PfCf dP ^ _ 
kf dt 
P(r,0) =0, 
P{oo, t) = 0, 
, dP 

-2'Kkf — 
or 


—( 

dP 

dr V 


r > 


t>0 

— 

1. 


r > rp, t > 0 




( 6 ) 


{pf and Cf are mass density and heat capacity of the soil, respectively). 
It can be shown as in van Everdingen and Hurst 0 that p{t) = 0{Vt) 
and p'{t) = 0{l/\/t), as t —> 0. Therefore, p'{t) is in L[q^([0,oo)) and 
we can make partial integration in ©. Taking natural assumption that 
qf = 0 at t = 0 (that is T^. = at t = 0) and using (|1J) we obtain a 





4 


Volterra integral equation for qf. 

bf{Te{z,t)-T^{z)) = qf{z,t)+bf [ p {t - T)qf{z,T) dr. (7) 

Jo 

This equation has the resolvent r G Lj'Q^([0, oo)) and it can be solved by 
formula (see Gripenberg, Louden and Staffanson |3]) 

= bf[Te-T^-r-k (Te - T^)], (8) 


where we have introduced convolution operator 

{r*(j)){t)= f r{t — t)(I){t) dr. 

Jo 

By use of (jH| and 0 we can eliminate formation temperature from © 
which is then transformed to 

8T BT B‘^T 

ae[^ - Ve-^] + be{T, - Ti) + bf{T, -r^T,) = + F, (9) 

where F = bf(Tz — r -kT^) is a smooth known function. We see that 
equations du, m and © represent parabolic system perturbed by the 
operator M given by 


Mu{z,t) = f r{t — t)u(z,t) dr (10) 

Jo 

The problem is to solve the system o, m and © with suitable 
boundary and initial conditions. We assume given the temperatures of 
entering water at z = 0 and oil at 2 : = L. At the bottom of inner and 
outer annulus we have equality of water temperatures and continuity of 
total thermal flux. Therefore, we take 


BT BT 

-^{0,t) = 0, To{L,t)=T^, Ti{0,t)=Tj^, 
BT BT 

Ti{L,t)=Te{L,t), D,-^{L,t) + De-^{L,t) = 0, 


( 11 ) 

( 12 ) 


for all t > 0, where T^ = Tz{L) and Tf are given. The initial conditions 
are 


To{z,Q) =Te{z,^) =Tz{z), Ti{z,^) =Tl{z), (13) 

where function satisfies compatibility conditions T^(0) = Tf, Tl{L) = 
Tz{L) and it is close to geothermal temperature T^. All functions in¬ 
volved are supposed to be smooth. 
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2. Variational problem 

We consider variational formulation of the problem o, m and © 
with boundary and initial conditions CUi-dlSI). Without lose of gener¬ 
ality we can consider homogeneous boundary conditions = Tf = 0. 
We introduce Hilbert space 

H = {(</.,, Hi(0,L)3 :(/)o(L) = 0, = 0, UL) = UL)] 

with the norm || • || inherited from H1(0, L)^ and bilinear forms A, B and 
C over V xV defined as follows: for T = (To, Ti,Te), $ = {4>o, 4>i: 4>e) we 
set 


>t(T, $) — Ao{To, (fo) + Ai{Ti, (pi) + Ae{To, cfe) + B^T, $) 
Ao[To, (po) = j - aoVo-^f)o) dz 


Ai 


° dz dz 

dTi dePi 


dz 

dTi 


i{Ti,(pi) = J 

A(T^^ [\r.dT,d(P, 9T 

Ae{Te,(pe)=J [0^ —— - aeVe^(pe) dz, 


H(T, $) = bo r{To - Ti){<Po - <pi) dz + he ^iTi - Te){(Pi - <Pe) dz 
Jo Jo 

+ bf TeCpedz 

Jo 

CiT,^) = -bf [ {r^Te)(Pedz. 

Jo 

Duality between V' and V will be given by the formula 

(F, (Iq{Fo, (po) T Oii{Fi: (pi) T Q.e(-^e5 4^e) 

where F e W is of the form F = {Fo, Fi, Fg), Fo,Fi,Fe € {H^{0,L)y, 
and brackets at the right hand side signify duality between (Hi(0,L))' 
and H1(0,L). We set H = L^(0,L)^, with usual norm denoted by 
I • I, and by identifying H with its dual we have V C H C V', with 
dense and continuous injections. Furthermore, by WiV^V) we denote 
the space of all functions from L^(0, tmax;F) with time derivative in 
-h^(0, tmax! y')- It is Well known that W (F, V) is continuously embedded 

in C'([0,tmax];H). 

With this notations we can reformulate the problem O, ©, ©, 
(fTT |) -(P )) in the following variational problem: find T G W{y,V') such 
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that T(0) = £ H and for a.e. t £ (0,tmax) 

(T',$)+^(T,$)+C(T,$) = (F,$), y^£V. (14) 

The linear form on the right hand side is given by 

(F,$)= [ F<Pedz. 

Jo 

and it is obviously continuous. 

It is easy to see that .A(-, •) is continuous bilinear form on V which 
satisfy 

^(T,T)+7|T|2 > a||T||2, VT e F, 

with some constants a, 7 > 0. Bilinear form C{-,-) comes from per¬ 
turbation operator M. It is not difficult to see that for any function 
W. ( 0 , tmax ) —> L^( 0 , L) it holds 

\\Mu{t)\\L2(^o^L) < k(i-'r)||k(T)||L 2 ( 0 ,L)dr)^/^ (15) 

where JC{t) = J^\r{T)\dT. From here it follows that M is linear and 
continuous operator from L°°( 0 , tmax; .^^( 0 , L)) to itself, and it has the 
following continuity property: if Un-,u £ L°°( 0 ,tmax; L 2 ( 0 ,L)) are such 
that 

Un{t) u{t) in L^( 0 ,L) for a.e. t £ ( 0 ,tmax) 

then 

Mun{t) —> Mu{t) in L^(0,L) for a.e. t £ (0,tmax)- 

Furthermore, it is easy to see that M is an operator of local type, as 
dehned in Artola PP, and therefore we can apply Theorem 1 from ^ 
and conclude: 

Theorem 1 Variational problem has a unique solution T £ W(V, V') 
for any £ H and F £ L^(0, tmax; y')- 


3. Numerical approximation 


In this section we discuss numerical approximation by finite difference 
method of the problem O, ©, ®, (EB-CSl). Instead of using equation 
0 we find more convenient to apply finite difference method to the 
equations O, 0, 0 and to discretize directly equations and 0 . 
We avoid numerical resolution of problem 0 by the use of Hasan and 
Kabir approximation: 


P{t) = Pn{ 
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where 

^ \i(0.80907 + log(s))(l + ^) fors>1.5. 

Furthermore, in our problem constants Dq, Di and Dg are very small 
and it is natural to consider hyperbolic system = Di = Df, = ^) 
instead of parabolic one. Due to limited space we will not enter here into 
discussion of existence theory for hyperbolic system. We just note that 
any difference scheme adapted to hyperbolic version of the system o- 
m will produce certain amount of numerical dispersion that will cover 
thermal diffusion in equations m-m, at least for reasonable mesh sizes. 
Therefore we chose to neglect thermal diffusion and consequently to drop 
superfluous Neumann boundary conditions for oil and water in outer 
annulus. This will generally change the solution just in corresponding 
boundary layers. 

We apply explicit finite difference scheme of first order with convective 
terms treated by upwinding. In all the experiments we have used a 
uniform grid in space and time. The spatial step h and time step r are 
related by the fixed positive number A through relation \ = t jh. 

In discretization of integral equation © we use composite trapezoidal 
rule which gives the following procedure for calculation of formation 
temperature at t = nr and 2 : = ih'. 


n—1 


rpn _ 

-^FA - 


i + Pi 


Tz,i + - T!^,i)(Pn+l-k - Pn-l-k) + PlTv,i 


k=l 


As a consequence of the convolution in formula © we see that the 
solution on next time level includes the solutions on all previous time 
levels. 

It can be shown that that described explicit scheme is TVB (total 
variation bounded) and L°°-stable if the following CFL condition is sat¬ 
isfied: ^ 

“ max{vo,Vi,Ve} + Ch' 

where C > 0 is certain constant that can be calculated from the coeffi¬ 
cients in ©-©• 

We now proceed with some numerical results. To evaluate the merits 
of one flow arrangement over another (counter flow and parallel flow), 
some conditions must be equal. The interval of time during which the 
water is cooled is not equal to the interval of time during which the water 
is heated. The sum of these time intervals we call circulating period or 




cycle. Both method can now be compared using the same circulating 
period. 

Results of our simulations after four cycles are presented in the figure 
Fig. 2. Counter flow heat exchange temperature calculations are shown 
on the left figure. The tubing temperature is almost always less then 
inner and greater than outer annulus temperature. 

Parallel flow heat exchange temperature calculations are shown on 
the right figure. The tubing temperature lies between the inner annulus 
temperature and formation temperature. In this case oil temperature is 
lower than any annulus temperature. Besides, formation temperature is 
higher than in the previous case. 

Tubing temperature as well as outer annulus temperature reach very 
soon almost constant level. The important thing to note with respect to 
the bottom-hole fluid temperature is that this temperature continually 
changes with time. A steady-state condition is never attained. Hence 
the stabilization of both outlet temperatures does not mean that all of 
the temperatures in the circulating system are constant. 

Under the same conditions we found that in parallel-flow arrangement 
temperature drop is smaller. Therefore, we may conclude that parallel 
flow seems to be better. 




Figure 2. Temperature calculation, left for counter flow heat exchange and right for 

parallel flow heat exchange. 

Legend: !=□ inner annulus, 2=0 outer annulus, 3=0 tubing, 4=V earth 

To conclude we point out that the linear model presented in this article 
has simplicity as its main advantage. It is not difficult to implement it 
in a computer code and it gives certain initial estimate of heat exchange 
in the system. Yet, important physical processes, such as dissipation 
due to friction and variations of viscosities and mass densities with the 





Heating of oil well by hot water circulation 


9 


temperature, are not taken into account. They lead to nonlinear model 
that will be considered in our forthcoming publication. 
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